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DISCLOSING THE RADIO LOUDNESS DISTRIBUTION DICHOTOMY IN QUASARS: AN UNBIASED 
MONTE CARLO APPROACH APPLIED TO THE SDSS-FIRST QUASAR SAMPLE 

M. Balokovic 1,2 , V. Smolcic 2,3,4 , Z. Ivezic 5 , G. Zamorani 6 , E. Schinnerer 7 , B. Kelly 8 

ABSTRACT 

We investigate the dichotomy in the radio loudness distribution of quasars by modelling their radio 
emission and various selection effects using a Monte Carlo approach. The existence of two physically 
distinct quasar populations, the radio-loud and radio-quiet quasars, is controversial and over the last 
decade a bimodal distribution of radio loudness of quasars has been both affirmed and disputed. 
We model the quasar radio luminosity distribution with simple unimodal and bimodal distribution 
functions. The resulting simulated samples are compared to a fiducial sample of 8,300 quasars drawn 
from the SDSS DR7 Quasar Catalog and combined with radio observations from the FIRST survey. 
Our results indicate that the SDSS-FIRST sample is best described by a radio loudness distribution 
which consists of two components, with 12 ± 1% of sources in the radio-loud component. On the other 
hand, the evidence for a local minimum in the loudness distribution (bimodality) is not strong and we 
find that previous claims for its existence were probably affected by the incompleteness of the FIRST 
survey close to its faint limit. We also investigate the redshift and luminosity dependence of the radio 
loudness distribution and find tentative evidence that at high redshift radio-loud quasars were rarer, 
on average louder, and exhibited a smaller range in radio loudness. In agreement with other recent 
work, we conclude that the SDSS-FIRST sample strongly suggests that the radio loudness distribution 
of quasars is not a universal function, and that more complex models than presented here are needed 
to fully explain available observations. 

Subject headings: Galaxies: surveys - Cosmology: observations - Radio continuum: galaxies 



1. INTRODUCTION 

An important property of Type 1 AGNs (i.e. broad line 
quasi stellar objects; QSOs) is the existence of radio-loud 
(RL) and radio-quiet (RQ) populations. One of the more 
controversial topics in QSO studies is whether these RL 
and RQ quasars^ form two physically distinct popula- 
tions of objects. On the one hand, both types of quasars 
are likely powered by similar physical me chanisms (e.g. 
Barth el 1989; Urry & Padovani 1995; iShankar et all 
2010!) , and their radio loudness has been shown to anti- 
correlate with accretion rates (in Eddington luminosity 
units) onto t heir c entral supermassive black holes (e.g. 
iSikora et al.l [2007). On the other hand, it has been 
demonstrated that, relative to RQ quasars, RL quasars 
are likely to reside in more ma ssive host galaxies (e.g. 
Peacock, Miller, & Longair 1986: ISikora et al.ll2007l ). and 
harbor more massive central black holes (e.g. Laor 2000; 
Lacy et al. 2001; McLure & Jarvis 2004). However, 
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mainly due to two problems, i) severe selection biases, 
such as incompleteness in observed QSOs samples, and 
ii) the overwhelmingly high fraction of RQ quasars, it 
is still unclear whether RL and RQ quasars form two 
distinct populations of objects, or a c ontinuous sequence 
(as suggested by e.g. lLacv et al~ll2001[ ). Here we focus on 
this problem. 

Whether RL and RQ quasars form two distinct pop- 
ulations can be studied by investigating the relation 
between their radio and optical emissions. This illu- 
minates the relative importance of the likely dominant 
sources of electromagnetic radiation in these two wave- 
length windo ws, namely the relativist ic jet and the ac- 
cretion disk. iStrittmatter et all (|1980l ) first pointed out 
that the radio-to-optical flux density ratio for optically 
selected QSOs appears bimodal (the so-called 'quasar 
radio dichotomy'), which suggests that QSOs can be 
divided into two distinct, RL and RQ, populations of 
objects. Other a uthors found addi t ional evidence for 
a dichotomy (e.g . iKellermann et al.l 119891 : iMiller et al. | 
| 1990tllvezic et al.l 120021 : iWhite et aUl2007t iZamfir et al.l 
I 2008D. However several studies (e.g. IWhite et al.l l2000t 
iCirasuolo et al.ll2003bt lLacv et al.ll2001l ) disputed its ex- 
istence. A bimodal distribution in quasar radio loudness 
would point to distinct physical properties of radio-loud 
and -quiet QSOs, such as a d ifferent physical orig in of ra- 
dio emission (jet vs. corona; lLaor fc B char 2008), differ- 
ent black hole masses (e.g.lMcLure fc Jarv is 2004) , accre- 
tion rates ( e.g. ISikora et al.ll2007tpamiltonll2010l) and/or 
spins Ce.g. iBlandford fc Znajekl 119771 : iBlandforcfl fl990l : 
iGarofalo et al.l I2010D. as well as host galaxy properties 
(e.g. ISikora et al.l 120071: lLagos et al1l2009t Kimball et al.l 
120111 ). The existence of two distinct quasar populations 
may then be linked to hierarchical structure growth in a 
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ACDM universe in which dark ma tter halos and galaxy 
mergers play an important role (e.g.lHuehes fc Blan dford 
2003 : iLagos et all 1 2009; Ga rofalo et all 120101: Hamilton! 
2O10t iFanidakis et alJl2010t l. 

The radio-to-optical flux density ratio as a measure of 
radio loudness was initially proposed by Schmidt (1970). 
A division between RL and RQ quasars at a ratio of 
~10 was found by Kellermann et al. (1986) using rest- 
frame 5 GHz (6 cm) radio and B-band optical flux 
densities. An alternative radio loudness definitional is 
purely based on radio lu minosity, e.g. L 14 GHz > 10 25 W 
Hz" 1 (jMiller et al.lll990D . Studies of large QSO samples 
have yielded that the most RL QSOs are ~ 10 3 times 
more luminous in radio than in optical (e.g. Ilvezic et al.l 
2002). On the other hand, much deeper radio stud- 
ies of smaller samples showed that RQ QSOs typically 
have a factor of ~ 10 weaker radio than optic al emission 
/e.g. IKellermann et alJlT98flh. Using SPSS (lYork et al.l 
2000h and F I RST (IBecker. White fc Helfandlll995l) data. 
Ivezic et al.l (|2004bl ) have shown the existence of a peak 
at the high end of the radio-to-optical flux ratio distri- 
bution. However, because of the detection limits of both 
surveys, they could not detect a significant number of 
RQ objects at radio wavelengths, and hence argued the 
existence of a secondary peak at the low end (and thus 
bimodality) based on two arguments: i) the majority of 
SDSS quasars (~ 90%) were not detected in the FIRST 
survey, thus they should lie at the low end of the radio 
loudness distribution, and ii) deep radio studies showed 
that typical RQ QSOs have a factor of ~ 10 4 weaker 
radio emission than typical RL QSOs, thus a secondary 
peak must exist. 

One of the most recent results on the quasar radio di- 
chotomy comes from a stacking analysis of FIRST 20 cm 
snapshot images at the posit ions of ~ 40 000 quasars 
from the SDSS DR3 catalog (IWhite et al.ll200l . They 
showed that a shallow minimum exists between the RL 
and RQ parts of the radio loudness distribution; how- 
ever, it was stressed that optical selection effects proba- 
bly dominate the observed distribution of radio loudness. 
Specifically, the authors discuss the difficulty of identifi- 
cation of SDSS quasars at 2.4 < z < 3, w here their colors 
are v ery similar to stellar colors (see e.g.. lRichards et aTl 
2002) , and support their claim by noting the strengthen- 
ing of the bimodality for a subsample of quasars in that 
r edshift range. 

ICirasuolo et al.l (|2003bl ) have used Monte Carlo simu- 
lations of the quasar population in order to model the 
intrinsic radio loudness distribution. They have com- 
pared their simulated samples to three quasar samples 
(2dF QRS, LBQS and PBQS; see ICirasuolo"!^ l2l)03bl 
for detail s) with FIRST su r vey da ta and radio observa- 
tions by Kell ermann et all (|1989l ). The three samples 
probe a wide range of radio loudness, but contain fewer 
than 200 radio-detected quasars in total. The authors 
reported that the best-fit radio-to-optical flux ratio dis- 
tribution is a double-Gaussian function with ~97% of 
the quasars in the RQ, and ~3% in the RL component. 
However, they conclude that there is no minimum be- 
tw een the two Gaussia ns, thus no bimodality exists (but 
see Ilvezic et al.ll2004bl for a different interpretation). 

10 For more details on definitions of radio loudness the reader is 
referred to Appendix C in Ivezic et al. (2002). 



The main problems causing the discrepancies in litera- 
ture regarding the existence or absence of a quasar radio 
loudness dichotomy lie either in the ambiguity of quasar 
selection, low number statistics, or severe selection biases 
linked to flux-limited samples (because radio loudness is 
defined as a radio-to-optical flux ratio, i.e. a ratio of two 
quantities drawn from flux limited samples) . To properly 
address all of these biases, in the work presented here 
w e utilize a Monte Carlo based approach similar to that 
in ICirasuolo et al.l l2003bl but more robust and using a 
much larger sample of observed quasars (the largest QSO 
database currently available is the SDSS quasar catalog). 

The outline of the paper is as follows. In Section [5] 
we present the data and our sample. In Section [3] we 
outline our methodology and the simulation algorithms. 
We present our results in Section 01 and discuss them in 
Section EJ We summarize our conclusions in Section |6l 
Throughout the paper we report magnitudes in the AB 
system, for the optical as well as for the radio. We use 
standard cosmology (H = 70 km s _1 Mpc -1 , Qm = 0.3, 
Oa = 0.7) and quasar continuum spectrum defined as 
/„ cc v a , where v denotes frequency, /„ is flux density 
and a is the spectral index. 

2. DATA AND SAMPLES 
2.1. Choice of the Source Catalog 

For a robust study of the radio loudness of quasars it 
is necessary to use a large, statistically significant radio- 
optical sample of quasars that covers the broadest possi- 
ble range in radio loudness. As quasars are rare objects 
(e.g. the space density of SDSS quasars with i < 19 is 
~11 deg." 2 ) a statistical sample can be assembled via 
large-area observations. On the other hand, achieving 
radio sensitivities over these fields deep enough to probe 
far into the radio-quiet regime with the current genera- 
tion of radio interferometers is challenging, and unfeasi- 
ble over fields larger than a few degrees. 

Because of the scaling of both radio and optical depth 
limits with surveyed areas, state-of-the-art surveys, such 
as e.g. SPSS-FIRST (~9380 deg. 2 , % < 22.5, > 
1 mJ y; iSchneider et al.ll2010t IBecker. White fc Helfandl 
Il995lh Stripe 82 (92 de g. 2 , g < 24.5, S 1AG Hz > 
52 uJv: IHodge eTalll20Tl . COSMOS (2 deg. 2 , i < 26.5, 
Si 4GHz > 50 /xJy; Scoville et al. 2007; Schinnerer et al. 
2007), and VVPS (1 deg. 2 , I < 24, ft 4G Hz > 80 /zJy; 
iLe Fevre et 111120031 : iBondi et al.ll200l cover a compara- 
ble radio loudness range, while smaller area surveys suffer 
from small number statistics in their optical quasar sam- 
ples (see e.g. Smolcic et al. 2008). Hence, for the analysis 
presented here we have utilized the largest available sam- 
ple of quasars with radio coverage, drawn from the SDSS 
and FIRST sky surveys. 

For the reasons outlined above, we use the SPSS 
Seven th Pata Release Quasar Catalog (Schneide r et al.l 
2010). It contains 105,783 spectroscopically confirmed 
quasars taken from ^9,380 deg. 2 of the sky. The cat- 
alog consists of objects with spectroscopy that 1) have 
reliable redshifts, 2) have at least one emission line with 
FWHM greater than 1000 km s" 1 or interesting/complex 
absorption features, and 3) are more luminous than rest- 
frame Mi = —22.0 (in our cosmological model, assuming 
a power-law continuum with a = —0.5). The selection is 
based on sources' i magnitude and position in the multi- 
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dimension al color-space based o n five SDSS photometric 
bands (see iRichards et aLl[2002l for details). The quasars 
in the sample have 15.0 < i < 21.2, but the majority of 
quasars are brighter than i w 19 because the flux limit 
for the main spectroscopic sample is i < 19.1. In addi- 
tion to the main spectroscopic sample, the catalog con- 
tains serendipitously identified quasars and sources se- 
lected for their proximity (within 2") to a 1.4 GHz radio 
source drawn from the FIRST survey. The rest of the 
quasars in the sample are matched to the FIRST catalog 
to within a 2" radius. The radio flux densities at 20 cm 
(1.4 GHz) are given in "AB radio magnitudes": 

where / rad is the radio flux density at 20 cm in Jy. 

FIRST is a radio survey at 1.4 GHz/20 cm, conducted 
with the Karl G. Jansky Very Large Array' (V LA) in 
B-configuration (|Becker. White fc Helfandlll995l ). It has 
mapped approximately 10,000 deg. 2 of the North Galac- 
tic Cap with a beam size of 5.4" and a typical RMS sen- 
sitivity of 0.15 mjy/beam. Flux density values used here 
are integrated over the two-dimensional Gaussians fitted 
to each source. Due to the lack of very short spacings in 
the VLA B-configuration and the nature of the Gaussian- 
fitting source detection algorithm, fluxes for extended 
objects larger than about 10" are likely underestimated 
and split into multiple components in the FIRST source 
catalog. We estimate that this does not significantly af- 
fect our sample as most of the radio-detected quasars 
are expected to be unresolved at the angular resolution 
of the FIRST survey. Multiple component sources are 
rare and in most cases radio-loud, which puts them high 
above the interesting transition region between RQ and 
RL regimes (e.g. Jiang et al. 2007). For a discussion 
of the distribution of integrated and peak flux densities 
and the extended source bias of FIRST-detected SPS S 
quasars we refer the reader to Kimball fc Ivezic I (|2008[ ). 
Nonetheless, we do take this effect into account statis- 
tically, using the FIRST survey completeness correction 
derived specifically for SDSS quasars (see Figure 1 in 
Jiang et al. 2007, and references therein), taking into 
account both the source size and flux distribution. 

2.2. Our Main Sample 

Following the SDSS DR7 Quasar Sample notation, 
hereafter we adopt the apparent radio magnitude (t) 
for radio flux densities. The catalog lists 8,630 quasars 
with a radio detection within 2" from the optical 
source position. For 'radio luminosity' here we use the 
1.4 GHz/20 cm absolute radio magnitude (M t ), derived 
from the apparent radio magnitude using a K-correction 
of the form —2.5(1 + a) log(l + z) (e.g. Richards et al. 
2006), assuming a = —0.5 for each quasar (e.g. Kimball 
& Ivezic 2008). For the 'optical luminosity' we use the 
emission-line-corrected absolute magnitude in the rest- 
frame SDSS i band (A e// = 7471 A). The reason for 
this choice is that the SDSS i-band was used to con- 
struct the flux-limited DR7 quasar sample and that i 
magnitudes are available even for quasars at the highest 
redshift. Prior to computing the absolute i magnitudes 
for the quasar continua, we have corrected the cataloged 
apparent magnitudes for galactic extinction using the 



extinction maps of Schlcgcl, Finkbci ner fc Davisl (|1998|1 . 
To get continuum magnitudes, we subtract the contri- 
bution of emission lines t o the i band using a function 
derived by IRichards et al.l ([2006). Using the mean SDSS 
quasar spectrum, they have calculated the contributions 
of major quasar emission lines (above the power-law con- 
tinuum) as a function of redshift, up to z w 5. After 
the subtraction of this contribution to the «-band appar- 
ent magnitudes, we use the canonical K-correction with 
a = —0.5. Given the high uncertainty of this correction 
for z > 5, hereafter we exclude 56 quasars with z > 5 
from our sample. For d etails of the calculatio n see Sec- 
tion 5 and Figure 15 in Richards et al. (2006). In order 
to get a 'cleaner' optical selection, we exclude the quasars 
selected only on the basis of their proximity to a FIRST 
source. After the exclusion of those quasars, the z > 5 
quasars and those with uncertain photometry, we are left 
with 8,307 quasars in our main radio-optical subsample. 
Our main optical sample consists of 98,544 quasars. 

3. THE METHODOLOGY FOR MONTE CARLO 
SIMULATIONS 

If one had a complete volume-limited quasar sam- 
ple with both rest-frame optical and radio luminosities, 
studies of the radio loudness distribution would be sim- 
ple. For example, if radio and optical luminosities were 
linearly correlated (regardless of the physics of such a 
setup) , the radio-to-optical luminosity ratio distribution 
would be very narrow and the radio versus optical lu- 
minosity diagram would show a straight narrow band. 
Another possible situation could be a broad distribu- 
tion of the radio-to-optical luminosity ratio, with per- 
haps evidence for its two-component nature (dichotomy), 
with or without evidence for a local minimum (bimodal- 
ity). In this case, the radio versus optical luminosity 
diagram would show a broadly dispersed line or two pos- 
sibly overlapping lines, one for the RQ and one for the RL 
quasars (see the bottom two panels of Figure [T] for an il- 
lustration) . If the sampled redshift range was sufficiently 
broad, one could even perform such studies for subsam- 
ples selected from narrow redshift slices, and search for 
evidence of evolution in the inferred properties with cos- 
mic time. 

Unfortunately, such a complete sample does not exist. 
The main difficulty when working with real samples is 
that both optical and radio data are truncated at finite 
flux levels, and this strong selection effect must be prop- 
erly accounted for. There are additional effects, such as 
K-corrections and optical variability, which can be han- 
dled to some extent, as discussed further below. In order 
to utilize existing samples, various statistical methods 
have been used in published work to account for the ob- 
servational effects, as briefly reviewed in Section 1. 

Here we use Monte Carlo simulations, with the follow- 
ing main logical steps: 

1. The main sample is defined at optical wavelengths, 
with various selection effects, and the main observ- 
able is the optical apparent magnitude. The sam- 
ple also contains redshifts for all quasars and radio 
magnitudes for a subset of radio-detected quasars. 
The selection function for the sample is well known. 

2. Only a subset of real optical sources from the main 
sample are detected at radio wavelengths, compris- 
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ing the main radio-optical subsample. They have 
a certain distribution of radio magnitudes and, at 
least in principle, a different distribution of optical 
magnitudes than the main optical sample. The dis- 
tributions of radio and optical magnitudes for the 
radio-optical subsample are the main constraints 
on models for the relationship between the radio 
and optical luminosities. 

3. Given a quasar from the main optical sample with 
certain optical luminosity, a parametrized model 
generates its radio luminosity (as detailed below). 
Using its implied apparent radio magnitude and 
the radio selection function, this object is either re- 
tained or rejected from the simulated radio-optical 
sample. 

4. Starting with the main optical sample, the simu- 
lation produces a subset of quasars with simulated 
radio magnitudes - a simulated radio-optical sub- 
sample. The distribution of simulated radio magni- 
tudes for this subsample, as well as its correspond- 
ing distribution of optical magnitudes, is compared 
to observed radio and optical magnitude distribu- 
tions of the main radio-optical subsample, and uti- 
lized in a x 2 minimisation procedure to get the 
best-fit model parameters. 

Although the Monte Carlo sample utilizes observed 
optical magnitudes, its radio magnitudes are generated 
stochastically. Therefore, there is no object-to-object 
correspondence between the real radio-optical subsam- 
ple and the simulated radio-optical subsample. For a 
good model, we expect a statistical agreement for two 
main quantities: i) optical magnitude distribution, and 
ii) radio magnitude distribution between the observed 
and simulated radio-optical subsamples. For simulat- 
ing the radio luminosities we initially use four models 
independent of redshift and optical luminosity and the 
complete main sample presented in Section 12.21 How- 
ever, once we select the best model, we investigate the 
redshift evolution and dependence on optical luminosity 
for its parameters. The following subsection defines ra- 
dio loudness and describes the models considered in this 
work. The rest of this section describes the simulation al- 
gorithm (Section l3.2j) . the optimization strategy (Section 
13.31) and the evaluation of our method on purely artificial 
data (Section 13.41) . 

3.1. The Radio Luminosity Models 

We use four different models for the relationship be- 
tween the radio continuum luminosity at 1.4 GHz (L rad ) 
and optical continuum luminosity in the SDSS i band 
(£ opt ). Their relationship is parametrized using radio 
loudness (Rk) defined as: 

Rk = log (L tad )-K log (L opt ). (2) 

Here K is either or 1 and it selects one of the two dif- 
ferent families of models that correspond to two common 
definitions of radio loudness found in the literature. For 
K = the distribution of radio loudness simply equals 
the distribution of radio luminosities (thus RL quasars 
are defined as being m ore luminous than some threshold 
radio luminosity; e.g. iPeacock et al.l 119861 iMiller et al.l 



199C , also examined in llvezic et a l. 2002 and Jian g et all 
20071) . while for K — 1 radio loudness is defined as the 
logarithm of the radio-to-optical luminosity ratio (with 
RL qua sars having this ratio greater than some thresh- 
old: e.g.lKellermann et allll98allvezic et al.|[2002L l2004bt 
iCirasuolo et al.ir2003alj bl: I Jiang et al.ll2007l ). Each defini- 
tion is suitable for a specific model of the relationship 
between radio and optical luminosities. 

We have considered four simple models for the radio 
luminosity: two in which it is independent of the optical 
luminosity and two in which it is directly proportional 
to the optical luminosity. In the former case radio loud- 
ness is properly defined with K = 0, Rq = log(L rad ) 
and radio luminosity is modeled as L rad — 10 R °. In 
the latter case, the proper definition of radio loudness is 
K = 1, R\ = log (L rad /L opt ), hence the radio luminosity 
model is L rad — L opt x 10 Rl . In both cases, we examine 
single- and double-Gaussian distributions of radio loud- 
ness. The designations of the models in this paper and 
their basic descriptions (where PDF stands for probabil- 
ity density function) are: 
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Examples of all four models are plotted in Figure[TJ Mod- 
els Ml and M3 are two-parameter models, with the free 
parameters being x (the mean of the Gaussian) and a 
(the width of the Gaussian). Models M2 and M4 have 
five free parameters: two for each Gaussian (xi, a\, X2 
and 02) and an additional parameter /, which determines 
the ratio between the integrals of the two Gaussians. In 
all cases the overall normalization is determined auto- 
matically by the requirement that the number of radio- 
detected quasars in the simulated radio-optical samples 
matches the number in the observed SDSS-FIRST sam- 
ple. 

Initially, we use these four models to simulate radio- 
optical samples based on the complete main optical sam- 
ple of 98,544 SDSS quasars and match them to the radio- 
optical subsample. Later we search for redshift evolution 
and dependence of the radio loudness distribution shape 
on optical luminosity not by adjusting the models, but 
rather by separating the main sample into smaller sam- 
ples constrained in redshift and optical luminosity. Possi- 
ble trends can then be inferred by observing how best-fit 
model parameters change between different subsamples. 

3.2. The Simulation Algorithm 
3.2.1. The Main Optical Sample 

For all radio loudness models that we have considered, 
the simulated radio magnitude (i) was calculated for each 
optically detected quasar in the main optical sample pre- 
sented in Section l2~2l This sample contains 98,544 op- 
tically detected quasars with all of the selection biases 
introduced by the DR7 Quasar Catalog, except for the 
radio-only selection which was eliminated from the cata- 
log for this work. Since we compare our simulated radio- 
optical subsamples to the radio-optical subsample of that 
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Fig. 1. — Examples of the radio models (Ml-4) considered in 
this paper plotted in the Mi-Mt plane. The radio loudness distri- 
bution is a single Gaussian for models Ml and M3 and a double 
Gaussian for models M2 and M4. Note that, in general, the Gaus- 
sians forming the two-component radio loudness distribution may 
exhibit complete, partial or no overlap, regardless of the model. In 
models Ml and M2 the radio absolute magnitude (Mt) is indepen- 
dent of the optical absolute magnitude (Mi), while in models M3 
and M4 Mt is proportional to Mi for every quasar. For clarity axis 
labels are shown only for the bottom-left panel, but they are the 
same for all panels. The dotted lines show where Mt = Mi. 



catalog (i.e. to the sources also detected in the radio), this 
biased sample is indeed the most valid optical sample for 
our simulations. 

Another option for the optical sample would be a sim- 
ulated sample produced from an empirical luminosity 
function. However, such an approach suffers from two 
main problems: i) a simulated optical sample would be 
free of selection biases (assuming selection bias was prop- 
erly corrected in the construction of the luminosity func- 
tion) and therefore we would additionally need to simu- 
late the SDSS optical selection; and ii) observed optical 
counts are difficult to simulate with proper uncertainties 
included (e.g. from photometry uncertainty, conversion 
between photometric bands, K-corrections, optical vari- 
ability etc.) unless the same uncertainties have been con- 
sidered during the generation of the luminosity function 
(see e.g. lLa Franca fc Cristia ni f 997 for an approach that 
does involve substantial consideration of uncertainties). 

3.2.2. Assignment of Apparent Radio Magnitudes and 
Application of the Radio Selection Function 

In each simulation, we assign an apparent radio mag- 
nitude (t) to each quasar from the main optical sample. 
The apparent radio magnitude t is calculated from the 
absolute radio magnitude M t , which is modelled accord- 
ing to Equation © as 

M t = KM, - 2.5R K , (3) 



where A takes on values of or 1 and Rk is a ran- 
dom variable drawn from either a Gaussian or a double- 
Gaussian probability distribution, depending on the 
model (see Section l3~T]) . 

To convert from M t to t we use a K-correction of the 
form -2.5(1 + a)log(l + z) (e.g. Richards et al. 2006), 
assuming a — —0.5 for each quasar (e.g. Kimball & Ivezic 
2008), and including corresponding uncertainties. Thus, 

t = M t + DM(z) + A c , rad (z) + n(a t ), (4) 

where DM(z) is the distance modulus and A c , rad (z) is the 
K-correction for a — —0.5, both dependent on quasar's 
rcdshift. The last term is a normally distributed random 
variable accounting for uncertainties in K-corrections and 
photometry. Its mean is zero and the standard deviation 
(at) can be estimated by considering the scatter in the 
meas ured radio spectral indices (e.g. iKimball fc Ivezic 1 
2008). As the standard deviation in spectral indices 
is approximately 0.3, it follows from the definition of 
the K-correction that a t ~ 0.35 for z = 2, where the 
number density of quasars is maximal. This uncertainty 
completely dominates the small photometric uncertainty 
(typically 0.03-0.04 mag). 

For our A = models (MI and M2) Equation © 
simplifies to M t = — 2.5i?o, hence 

t = -2.5A + DM(z) + A c , rad (z) + n(a t = 0.35). (5) 

For our A = 1 models (M3 and M4) 

t = Mi-2.5Rx + DM(z) + K c , t!ld (z) +n(a t =0.35). (6) 

The absolute optical magnitude is 

Mi = i - DM(z) - A c , opt (z) + n(o-i) + e(a v ), (7) 

where DM(z) is the distance modulus, A copt (z) is the 
K-correction and the remaining two terms account for 
the scatter in K-correction, photometry (normal dis- 
tribution, n(crj)), and variability (exponential distribu- 
tion, e(a v )). The optical spectral index distribution 
with an average of —0 .5 and a standard deviation of 0.3 
(|Richards et al.ll2006[ ) yields <7j ~ 0.35, taking into ac- 
count that the optical photometric errors (typically 0.03 
magnitude) are negligible. Optical variability scatter 
has an empirically determined exponential distrib ution 
with zero mean and a v » 0.2 ([Ivezic et~aT1[2004al) . Fi- 
nally, inserting the expression for Mj into Equation ([6]), 
the distance moduli cancel out, as well as the mean K- 
corrections (because of our assumption of the same mean 
spectral index in the optical and radio bands), and the 
normally distributed uncertainties at and ai add up in 
quadrature, yielding a convenient expression for t in our 
A = 1 models: 

t = i- 2.5Ai + n(a = 0.5) + e(a v = 0.2). (8) 

Our simulation algorithm assigns radio magnitudes 
to optically detected quasars by drawing random num- 
bers from appropriate distributions specified by Equa- 
tions ©, for Ml and M2, and ©, for M3 and M4. 
Each simulated radio-optical sample (consisting of 98,544 
quasars with apparent optical and radio magnitudes) is 
then subjected to the radio selection function in order to 
produce a smaller subsample of 'radio-detected' quasars 
corresponding to the observed radio-optical sample from 
SDSS and FIRST. All quasars fainter than the FIRST 
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flux limit at 1 mJy are rejected. Some of the quasars 
above the FIRST flux limit are randomly rejected in or- 
der to mimic survey incompleteness. This is performed 
by randomly choosing and excluding a number of quasars 
in narrow f-magnitude bins, with the fraction of ex- 
cluded objects being determined by the completeness 
fu nction of the surv ey up to its flux limit (see Figure 1 
in lJiang et al.ll200"7j) . At the end of this procedure, each 
simulated radio-optical subsample is one particular real- 
ization of the real radio-optical subsample drawn from 
SDSS and FIRST. The simulated radio-optical subsam- 
ples after radio selection have to consist of approximately 
8,300 quasars to at least roughly match the observations, 
which sets the overall normalization for our models. 

3.2.3. Evaluation of Goodness of Fit for the Simulated 
Samples 

Each of the simulated radio-selected subsamples ("re- 
alizations" hereafter) is binned into 10 bins in t, i, z and 
R' = 0A(i — t) distributions. We chose to examine these 
four distributions to provide better constraints for our 
models; in principle, weaker constraints, always consis- 
tent with those derived by using all four distributions, 
can be can be obtained with any subset of these distri- 
butions that includes either t or R' . Note that the R' 
distribution is a proxi to the radio loudness distribution 
in the case of K = 1; otherwise it is just an additional 
constraint on the relation between t and i magnitudefEH- 
The bin values are added to a pool of realizations simu- 
lated with identical model parameters. The radio magni- 
tude assignment, the selection procedure and the binning 
need to be repeated a large number of times for a given 
model and a set of its parameters in order to properly 
account for the stochasticity of particular realizations. 
The aim of this Monte Carlo procedure is to calculate 
the mean distributions of t, i, z and R 1 for a given set 
of model parameters and derive their expected variance. 
We calculate the mean and the standard deviation of the 
number of simulated objects in each bin (SNj and (jj in 
the equation below) from a large number of realizations 
and compute the total x 2 with respect to the binned t, 
i, z and R' distributions of the real SDSS-FIRST data. 
The x 2 is defined as 



X 



3 = 1 V 



SNj 



(9) 



where the sum runs over all 40 bins, while RNj is the 
number of real quasars in a particular bin and SNj and 
Oj are the mean and the standard deviation for that bin, 
determined from the simulations as described above. 

Although the four binned distributions are not entirely 
independent of each other (in particular, R' is a linear 
combination of i and t), the tests on artificial data (see 
Section I3.4[> have shown that the number of degrees of 
freedom may be assumed as if the bins were indepen- 
dent. Therefore, the number of degrees of freedom equals 
N — p — 1, where N is the total number of bins used for 
evaluation and p is the number of free parameters of a 
particular model. Hereafter we refer to one evaluation 

11 Also note that this distribution is severely biased and cannot 
be directly used to infer the distribution of radio loudness; see 
Ivezic et al. (2002, 2004) for an explanation. 
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to a unique value 



(represented here as the mean \ 2 for N r = 3000) as the number 



An illustration of convergence of \ 
X 

of realizations per simulation (N r ) increases. While smaller N r re 
quires less computing time, such simulations have greater variance 
of x 2 between them. For each given number of realizations per 
simulation, the same model (M4 here, with all parameters fixed) 
was simulated 1000 times in order to find the mean (large black 
rectangles), the standard deviation (smaller black rectangles at the 
ends of black lines) and minimum/maximum values of x 2 (g ra y 
rectangles at the ends of gray dotted lines) for a given N r . The 
insets show distributions of x 2 values for three selected N r values 
indicated in the upper left corner of each inset. The dashed ver- 
tical lines in the insets show the mean (thick) and the standard 
deviation (thin lines). 

of x 2 from its distribution given a large number of real- 
izations, as a "simulation" and abbreviate the number of 
realizations per simulation with N r . 

It is a general property of Monte Carlo simulations 
to converge to an average result only after the experi- 
ment is repeated a large number of times. We test the 
convergence by examining how the means % 2 , standard 
deviations and minimum/maximum values for different 
simulations with same N r . This is illustrated in Figure [2] 
for model M4 (analogous results follow from any of the 
models considered here) with a fixed arbitrarily chosen 
set of its parameters and a range of different N r . Clearly, 
a larger number of realizations per simulation reduces 
the probability of getting outlying \ 2 values that signifi- 
cantly deviate from the mean x 2 f° r the simulated model, 
but at the expense of computing time. We have found it 
more time-effective to work with a relatively small num- 
ber of realizations in our optimization procedure and cor- 
rect for the scatter in the % 2 values just before deriving 
the marginal probability density functions for model pa- 
rameters (see Section I5T3")) . Due to computing time con- 
straints, we chose to calculate 1000 realizations per sim- 
ulation as a reasonable compromise between accuracy 
and computing time. The calibration shown in Figure 
[2] shows that for N r = 1000 one should expect varia- 
tions with standard deviation of ~2%. For N r = 3000 
the variation drops to 1.2%, which does not represent a 
significant improvement compared to N r = 1000. The 
computing time scales linearly with the number of re- 
alizations per simulation, so N r — 3000 would require 
three times more computing time than N r — 1000. 

3.3. The Optimization Strategy 
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Since there is no a priori indication that the param- 
eter space for our models is simple, e.g. that there is 
a single minimum of the total x 2 , we began to search 
the parameter space with a random walk across a wide 
range of parameters. The algorithm we used to initially 
sample the parameter spaces with 2 to 5 dimensions, is 
a Metropolis algorithm (Metropolis et al. 1953). Sev- 
eral runs of the algorithm were used to identify the area 
around the global minimum of x 2 for each of the models. 
The program would then proceed to evaluate the region 
on a regular grid, so that marginal probability distribu- 
tions could be estimated. We have used three levels of 
refinement of the grid in order to be able to sample well 
the narrow region around the global minimum of x 2 ■ 

Since variations displayed in Figure [2] are expected to 
occur in the evaluation process, we needed to take into 
account the fact that some simulations resulted in \ 2 val- 
ues as low as the minimum \ 2 value, although their mean 
X 2 would be much larger. For all parameter space points 
for which x 2 wa s found to be less than 5 a % 2 (recall that 
this is ^2% for N r = 1000) above the lowest value found, 
the optimization would proceed to map the surrounding 
parameter space in more detail (typically, with a factor 
of 3 more resolution in all dimensions). 

For computation of the marginal probability density 
functions for model parameters, we assigned them the 
probability normally associated with x 2 values: 



50 



p(x) oc e 



(10) 



where a? is a vector (a set of coordinates) in the param- 
eter space of some model. The marginal probability of 
a certain parameter value is then the sum of the proba- 
bilities over all other parameter space dimensions. The 
probability density functions were normalized a poste- 
riori so that the sum of probabilities of each bin in the 
parameter values equals unity. To cope with the problem 
of the internal scatter in x 2 values (see Section l3.2.3|) . we 
used a special procedure to compute realistic probabil- 
ity density functions which take into account that the 
computed x 2 may vary according to Figure [2] First, we 
estimate the standard deviation of x 2 from the calcula- 
tions illustrated in Figure[2]and described in the caption; 
e.g. the fractional standard deviation for N r = 1000 is 
~2%. Then, we calculate the corrected marginal proba- 
bility density functions by adding a normal random vari- 
ate with standard deviation equal to 2% of the lowest x 2 
value to the previously computed x 2 values, calculating 
the probability distributions in each case and repeating 
this procedure many times. This Monte Carlo method re- 
quires ^1000 repetitions for the resulting mean marginal 
probability density function to converge. We derive the 
asymmetric error bars on parameter values as intervals 
containing ±34% of the total probability around the me- 
dian of the marginal probability density function for each 
parameter. 

3.4. Verification of the Approach on Artificial Datasets 

In order to evaluate our method and better understand 
the possible uncertainties and biases it involves, we first 
simulated artificial datasets resembling real data in all 
aspects, differing only in the fact that their radio loud- 
ness distributions are exactly known. Treating them the 
exact same way as the observational data, we have re- 



20 optimizations with 5 parameters 
(total of 100 results) 



fit, (7=1.069 

ideal, a=1.0 




^derived Xreal /^derived 

Fig. 3. — The distribution of parameters derived from our 
parameter optimization procedure (Xj er j vec i) relative to the in- 
put parameters of our artificial datasets (X rea j), normalized to 
the derived standard deviation of the reconstructed parameters 
("derived)' '^' le data is displayed for 20 optimizations on artificial 
samples constructed with model M4, which has 5 free parameters. 
The vertical error bars represent Poissonian noise. The solid black 
line marks the Gaussian fitted to the data, which is a close match 
to the ideal Gaussian with a = 1 (dashed line), verifying that with 
our method it is possible to reconstruct the model parameters as 
expected from statistics. 

constructed their model parameter values and compared 
them to the input ones. 

For a set of artificial datasets constructed with the ra- 
dio models presented in Section [3.11 the parameters were 
reconstructed to within 3er uncertainties in all trials. The 
distribution of the reconstructed parameters, normalized 
to their respective standard deviations, is displayed in 
Figure [3l The lowest x 2 values reached in optimizations 
were x 2 ~ 35, yielding a reduced x 2 of ~1 if the prob- 
lem has the number of degrees of freedom as if the fitted 
distributions were independent (see Section l3.2.3[) . 

We have also tried fitting an M4 (K = 1) model to 
artificial datasets constructed with model M2 (K = 0). 
A good fit should not be possible in that case. The lowest 
X 2 values in most of these cases stayed above ~1000. 
The fitted parameters present in both models (e.g. x\, 
X2) would often be several standard deviations from their 
correct values. We consider this result to be an indication 
of what may be expected if none of the models considered 
in this work yield an adequate description of the observed 
data. Tests performed with other models are consistent 
with the above results. 

4. RESULTS 

4.1. Models Fitted to the Complete Main Sample 

We attempt to reproduce the observed SDSS-FIRST 
sample by performing simulations to the whole main 
sample described in Section 12.21 The results of these 
optimizations are given in Table [T] for all four models. 
Based on our analysis, the least likely models to fit the 
observations are Ml and M2, in which the radio luminos- 
ity of quasars is independent of their optical luminosity 
(K = 0). The minimum total % 2 values reached in our 
optimization procedure were ^4000 and ^2500 for mod- 
els Ml and M2, respectively. Each of the models has 37 
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TABLE l 

Best-fit median values and Ict uncertainties of model parameters for our four models, obtained from fits to the complete 
main sample of quasars. value of k is a propery of the models and parameters descri be t he gaussian peaks (xi and £2), 
widths ((71 and (t2), and their relative normalization (/). see section |3. ii for. details. 
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Fig. 4. — The marginal probability distributions for all 5 parameters of our model M4. The parameters x\ and a\ set the radio-quiet 
portion of the radio loudness distribution, / is the fraction of radio-loud quasars and X2 and <T2 set the radio-loud part. The distributions 
were computed through Mo nte Carlo simulations that compensate for the instability of the \ 2 computed using a low number of realizations 
per simulation (see Section l3.3l for details). 



degrees of freedom. For both models, the largest con- 
tribution to the total x 2 comes from the disagreement 
with the observed redshift distribution. They also fail 
to correctly reproduce the most populated bins in the 
radio and optical magnitude distributions. Model M2, 
however, fits the bright end of the radio magnitude dis- 
tribution considerably better and this is clearly reflected 
in the lower x 2 value. Values of x 2 generally fluctuate 
with a x 2 < 100 (see Figure [2] and Section [3. 2 .3(1 . making 
the difference of 1500 in total x 2 of the two models highly 
significant. 

Models M3 and M4, both involving a proportionality 
between the optical and the radio luminosity of quasars 
(K = 1), were found to provide better fits to the main 
sample compared to models Ml and M2. Their lowest 
X 2 values of ^2000 and ^1000, respectively, are signifi- 
cantly lower than for models Ml and M2. Both models 
can reproduce the observed redshift distribution much 
more correctly (especially the lowest redshift bin), but 
the overall match is not a statistically good one since 
both have only 34 degrees of freedom. The most signif- 
icant contributors to the total % 2 are the optical mag- 
nitude bins with the highest numbers of quasars, with 
the bright end of the radio magnitude distribution con- 
tributing considerably to the higher x 2 value for model 
M3. The non-optimal matching in the radio and optical 
magnitude distributions results in discrepancies in the 



R' = 0.4(i — t) distribution as well. 

4.2. The Best-fit Model for the Complete Main Sample 

Among the models considered in this paper, the lowest 
X 2 values for the complete main sample were achieved for 
model M4 (K — 1 and a double Gaussian radio loudness 
distribution). We plot the marginal probability distri- 
butions for the parameters of the best-fit model M4 in 
Figure [4] Note that the expected scatter in x 2 values 
between simulations (due to a limited number of real- 
izations) was taken into account prior to constructing 
the probability density distributions (see Figure [5] and 
Section 13.2.31 for details). The median values and la 
uncertainties derived from the marginal probability dis- 
tributions for all five parameters of the model are given 
in Table [TJ In Figure [5] we plot the t, i, z and B! dis- 
tributions for the observed SDSS-FIRST data (the main 
sample) and the median values obtained from a simula- 
tion (1000 individual realizations) performed with model 
M4 and the best-fit set of parameters. The error bars 
on each bin are the minimum and the maximum value 
occurring in that bin when parameters are shifted ran- 
domly within ±1(7 from their respective best-fit medians. 

In Figure [S] we show the radio loudness distribution 
for our best-fit model M4, as well as the associated la 
uncertainties. We compute the probability of bimodality 
of the best-fit radio loudness distribution by deriving the 
fraction of the M4 model parameters which result in a bi- 
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Fig. 5. — Distributions of t and i magnitudes, redshift and 
B! = 0A(i - t) for the observed SDSS-FIRST sample (thick 
gray histograms) and for a simulation with the best-fit model M4 
and parameter values set to the median values of their respective 
marginal probability distributions (thin black histograms; see Ta- 
ble [D. The error bars mark the highest and the lowest bin values 
achieved in simulations with parameters set la off their respective 
medians in all possible combinations (2 5 = 32 combinations). 
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4.3. Dependence of the Best-fit Model Parameters on 
Redshift and Optical Luminosity 

In order to investigate whether a certain parametriza- 
tion of dependence on redshift or optical luminosity could 
be added to our M4 model, we have divided our initial 
main sample into 4 bins in redshift and 4 bins in appar- 
ent magnitude and repeated the optimization procedure 
for each of them. Since most of these 16 subsamples span 
fairly narrow ranges in redshift and optical magnitude, 
the fits were done using only radio magnitude (t) and 
R' = 0.4(i — t) distributions and hence the number of de- 
grees of freedom for this case is 14 per subsample. The 
choice for the binning was such that every subsample 
contains ^6000 optical and ^500 radio quasars, which is 
still large enough to perform fits to a distribution divided 
into 10 bins. The binning is given in Tabled along with 
the results for each subsample. 

Each of the 16 subsamples was best fitted with model 
M4, although for the optically faintest subsamples model 
M3 was almost equally well fit. x 2 values for individual 
subsamples range from 10 to 32 for 14 degrees of free- 
dom and the lowest total \ 2 value summed over all 16 
subsamples is 291 for 224 degrees of freedom (reduced % 2 
is 1.3). This is much lower than the lowest x 2 value for 
the M4 model fitted to the complete main sample, which 
is 1053 for 34 degrees of freedom. This result further 
confirms that a simple M4 model, which represents the 
distribution of the radio-to-optical ratio as a universal 
function independent of redshift and/or optical luminos- 
ity, is not a satisfactory representation of the real quasar 
sample from SDSS and FIRST. 

The results are plotted in Figure [7] using median red- 
shift and absolute optical magnitudes of the subsamples, 
showing that there are no clear trends in parameters x\ 
and (7i and that possible trends exist in the remaining 
three parameters as a combination of dependence on red- 
shift and/or optical luminosity. Results for the optically 
brightest subsample of each redshift bin were used to plot 
the continuous change of the radio loudness distribution 
in Figure [H With the current results it is not possible 
to disentangle the two dependencies or to quantify them, 
but we do discuss the tentative trends in the following 
section. 



Fig. 6. — Plot of the radio loudness distribution for median values 
of the M4 model parameters (thick black line). The shaded area 
marks propagated ltr uncertainties on the parameters. The inset 
shows a zoom-in on the region where the radio-loud regime starts 
to dominate. 



modal distribution. We consider a distribution bimodal 
if there is at least a slight minimum between the RL and 
RQ regimes. We infer a probability of ^20% when we 
take the la confidence intervals of each parameter into 
account. Thus, we conclude that the likelihood of bi- 
modality in the radio loudness distribution of our best-fit 
model is relatively small. Note, however, that based on 
our simulation tests (Section l3.4[) the total x 2 value of our 
best-fit model (1000 for 34 degrees of freedom) indicates 
that it is not a statistically acceptable representation of 
the data. In order to find a better-fitting model, in the 
next section we test for possible dependence of the M4 
model parameters on redshift and optical luminosity. 



5. DISCUSSION 

5.1. Implications from the Parameters of the Best-fit 

Model 

Our best-fit model hints at the possibility that there 
might exist two distinct populations of quasars, assum- 
ing that the double-Gaussian parametrization is an ap- 
propriate one. Considering the fraction of radio loudness 
distributions that show bimodality within ler confidence 
intervals of our best-fit model parameters, the probabil- 
ity that the global radio loudness distribution is bimodal 
is ^20%. Therefore, we can neither confirm nor firmly 
exclude that the radio loudness distribution is bimodal. 
A more robust result, however, is that the radio loudness 
distribution of SDSS-FIRST quasars can be described 
much better with two Gaussians than with a single one. 
This would imply that RL quasars obey a different rela- 
tionship between radio and optical luminosity compared 
to RQ quasars. In this paper we parametrize our mod- 
els so that there is either no relationship between radio 
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TABLE 2 

The bin limits («min. *max. zmin, zmax), median optical magnitudes (imed and M^med) and median redshifts (zmed) for the 
16 subsamples used to examine changes in best-fit parameters of the m4 model with redshift and optical luminosity. 
Best-fit parameters of the M4 model and their associated la uncertainties are given for each of our 16 subsamples. 
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and optical luminosities (K = models, Ml and M2), 
or £ i.4ghz = ii-b.,„d x lO Rl (K = 1 models, M3 and M4) 
for both types of quasars. The latter parametrization, 
for which we find a significantly better fit than for the 
former one, implies that the constant terms of the re- 
lationship (parameters x\ and X2, locations of the two 
Guassian peaks in R\) and its scatter (a\ and 02, widths 
of the Gaussians in i?i) are different for the two classes 
of quasars. In general, the proportionality between the 
radio and optical luminosity might be different for each 
of the two quasar classes, but this type of a model was 
not considered in the work presented here. 

Investigating a possible dependence on redshift or op- 
tical luminosity in our SDSS-FIRST sample, we have 
found a statistically good fit when the full sample is di- 
vided into 16 subsamples, each fitted independently (re- 
duced x 2 =291/224=1.3, summed over the 16 subsam- 
ples). The parameters describing the RQ Gaussian (x\ 
and (7i) were found not to vary significantly, which is not 
unexpected since this Gaussian is constrained only by its 
fall-off towards the RL regime. A possible trend with the 
absolute optical magnitude is suggested for the param- 
eters describing the RL Gaussian (#2 and 02) and the 
radio- loud fraction (/). The optically bright quasars ap- 
pear to be better described by an RL Gaussian which is 
'louder' (larger X2), a smaller dispersion [02) and a lower 
radio- loud fraction (/) than the optically faint quasars. 
There is a slight possibility that the trends are biased 
by the redshift-luminosity correlation inherent in all flux- 
limited surveys, but our method was designed specifically 
to avoid that kind of bias. Our tests performed on artifi- 
cial datasets (see Section l3T4|) lend confidence that model 
parameters can be recovered correctly to within statisti- 
cal uncertainties under a variety of different conditions. 

Figure [8] shows a possibly significant change of the ra- 
dio loudness distribution shape implied from our results 
for quasars with i < 18.5 (the brightest subsample in 
each redshift bin). This apparent shift of the RL peak 
and the change in its width could be indicative that the 
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Fig. 7. — Each of the five panels in this figure shows the variation 
of best-fit parameters of the model M4 with absolute magnitude 
and redshift. Redshift and absolute magnitudes plotted are the 
medians of the subsamples (see Table [2). Colors and symbols mark 
different redshift bins. The dashed horizontal lines mark best-fit 
values for the complete main SDSS-FIRST sample and the dotted 
lines indicate la uncertainties. See Section 5.1 for a discussion of 
possible trends. 
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Fig. 8. — Radio loudness distribution as a function of redshift 
for flux-limited subsamblcs of quasars with i < 18.5. Note that 
bimodality seems to have been more prominent at high redshift 
and that it flattened out at low redshift. The same general trend is 
detectable in fainter subsamples, but the shape of the distribution 
does not change in a monotonous manner. 

radio loudness distribution becomes increasingly bimodal 
at high redshift. A monotonous change in the radio loud- 
ness distribution shape is clearly visible for our subsam- 
ples with i < 18.5, indicating that radio- loud quasars 
are rarer, 'louder' and less spread out in their loudness 
at high redshift. The same general trend is roughly de- 
tectable for subsamples with i > 18.5, but the change 
of the radio loudness distribution shape is not nearly as 
clear and continuous as for the i < 18.5 subsamples. We 
wish to emphasize that we understand selection effects 
for the i < 18.5 subsamplc much better than for the 
i > 18.5 subsample. First, the quasar targeting is com- 
plete only for i < 19 candidates. Second, the quasar 
variability will introduce an RMS scatter of 0.2-0.3 mag- 
nitudes and thus blur this supposedly sharp flux limit. 
Hence, our choice of i < 18.5 binning defines the faintest 
sample for which the simple SDSS optical selection func- 
tion is highly reliable. The lack of clear and continuous 
change in fainter samples may be reflecting the lower re- 
liability of the i > 18.5 objects in the SDSS quasar sam- 
ple, or indicate that the dependence of the radio loud- 
ness distribution on redshift or optical luminosity is not 
monotonous. We plan to investigate this in the future 
with more complex model and different statistics (e.g. 
maximum likelihood), so that more information may be 
extracted from the existing data. 

5.2. Comparison to Previous Findings 

Previous work on the issue of the bimodality in the ra- 
dio loudness distribution of quasars were either based on 
small samples of quasars, o r large but le ss reliable ones. 
For ex ample, iWhite et alj (|2000[ ) and iCirasuolo et al.l 
(|2003al fbh used spectroscopically confirmed samples of 
636, 141 and 195 radio - detecte d quasars, respectively, 
while llvezic et al.l (|2002L l2004b| ) used photometric sam- 
ples with ~4400 and ^10,000 radio-detected quasar can- 
didate sources. In comparison to those results, the re- 
sults presented here have considerably better statistics 
and reliability. Our result that a two-component model 
with a direct proportionality between radio and optical 
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Fig. 9. — A comparison of our r esults (thick soli d and d ashed 
lines) to the previ ous results by [Cirasuolo et al. (20033) an< 3 
llvezic et al.l l|2004bl) (dot ted and dot-da s hed gr ey lines, respec- 
tively). The results from [Cirasuolo et al. (2003b) were shifted by 
0.4 towards lower values to account for different optical bands, i 
versus B. The thick solid line is the intrinsic radio loudness dis- 
tribution (Ri) with parameters given in Table [T] for model M4. 
The thick dashed line shows the result one would obtain with our 
sample by not accounting for the incompleteness of FIRST survey 
near its flux limit. 

luminosities fits the obs erved data best is consistent with 
ICirasuolo et al.l ([2003b). who have used a method similar 
to ours on a much smaller heterogeneous sample. Using a 
stackin g analysis of t he FIR ST data to probe faint radio 
fluxes, IWhite et al.l (|2007t ) also found a strong depen- 
dence of the median radio luminosity of quasars on the 
optical luminosity. 

Our results can be compared to some of the previ- 
ous ones as plotted in Figure [S] The prominent min- 
imum separati ng the RQ and RL Ga ussia ns observed 
at i?i ~ 1.3 by ICirasuolo et all (|2003bD and llvezic et al.1 
(|2004b| ) is likely due to the incompleteness of the FIRST 
survey at its faint end. In the region around Ri ~ 1, 
FIRST is not more than ^70% complete for the faintest 
radio sources, but this was not taken into account in the 
earlier works. The plot also shows how our result would 
change in the case where we neglect FIRST incomplete- 
ness and assume it is 100% complete down to its flux 
limit - in that case we would have found a weak bimodal- 
ity with a minimum at R\ ~1. Its relative weakness in 
comparison to previous results can be understood as a 
difference in the fraction of optically faint quasars in the 
sample, which is lower in our case (45% here compared 
to ^75 in other cases), and which limits the influence of 
incompleteness to lower values of R\ . Note that the flux 
limit in the radio is the same for all three studies com- 
pared in Figure |H1 so optically fainter samples introduce 
incompleteness at a higher value of R\. Some of the 
difference with previous results may also be related to 
the uncertainties arrising from K-corrections and optical 
variability. 

A small degree of discrepancy exists in the fraction 
of RL quasars between our and previous results; / = 
(12 ± 1)% (i nferred here) compared to (8 ± 1)% from 
llvezic et all (|2004bD and (3 ± 2)% from ICirasuolo et al.1 
(|2003bD . We have defined the radio-loud quasar fraction 
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Fig. 10. — Dependence of the radio- loud fraction, defined with 
a threshold in radio loudness and calculated from a best-fit model, 
as a function of apparent i magnitude (grey symbols and lines). 
For comparison, the same result from Jiang et al. (2007) is plotted 
in black. M4 model parameters for each subsample were obtained 
independently and radio-loud fraction was calculated by integrat- 
ing the analytic form of the best-fit radio loudness distribution 
function (e.g. it is not the / parameter of model M4). 



as the ratio between the areas under the RL and RQ 
Gaussians, denoted in our models as parameter /. In 
other words, the fraction of RL quasars is the fraction of 
quasars whose radio loudness is determined by the RL 
Gaussian. Note that with this definition the range of ra- 
dio loudness between the RL and RQ peaks is shared by 
both types of quasars. iCirasuolo et al.1 ( l2003bl) have u sed 
that definition as well, but llvezic et alJ (|2002L l2004bl l de- 
fined radio- loud quasars as the ones having R\ > 1. This 
definition is more practical and a radio loudness may be 
unambiguously assigned to each individual quasar. With 
this definition, our RL fraction is (10 ± 1)%, consistent 
w ith the results of llvezic et al.l (|2004bD . 

Uiang et alJ (|2007t ) have found that the fraction of RL 
quasars depends both on optical luminosity and redshift, 
practically independent of the exact definition of an RL 
quasar. They have found that the fraction of RL quasars, 
when defined with a threshold in radio loudness is higher 
for optically bright quasars and at low redshift. Our re- 
sults tentatively confirm that such trends exist, albeit 
at low significance (see the lowest panel in Figure [7]). 
With our definition of the RL fraction (the fraction of 
quasars in the radio-loud Gaussian, /) the trend with 
respect to the optical luminosity seems to be just the op- 
posite - / tends to be lower for optically brighter quasars. 
Note, however, that even if the radio loudness distribu- 
tion shape changes (e.g. as in our i < 18.5 results shown 
in Figure [8j , the fraction of RL quasars defined with a 
threshold in radio loudness can remain constant or have 
an opposite trend than parameter / . In order to check for 
consistency with Uiang et al.l (|2007f) , we have performed 
additional fits to subsamples derived from the original 
main sample by dividing it into narrow bins in appar- 
ent optical magnitude. As displayed in Figure QUJ the 
fraction of RL quasars calculated from the best-fit M4 
models for each subsample, usi ng a threshold in radio 
loudness, matches the data from Uiang et al.l (|2007l) very 
well. 



In broad agreement with previous work by other au- 
thors, our current results imply that the radio loudness 
distribution did not have the same shape at all times in 
cosmic history. It might be possible to investigate this 
in the future using deep radio imaging of samples i n thin 
redshift slices. For example, iKimball et al.l (j2011h have 
used deep imaging with EVLA to probe a statistically 
complete sample of quasars in a narrow redshift bin at 
0.2 < z < 0.3. They found that the radio loudness distri- 
bution can be well explained with a superposition of two 
QSO populations (radio-loud, which is AGN-dominated 
and radio-quiet, dominated by star formation in the host 
galaxy), and their distribution of radio luminosity ap- 
pears similar to the best-fit two-component model pre- 
sented here. Since currently available survey data lacks 
either depth in the radio or significant sky coverage, the 
data from the newer generation of radio instrumentation 
(JVLA, SKA) will be essential in order to fully constrain 
the radio loudness distribution and ultimately, under- 
stand its physical origin. 

6. SUMMARY 

In this paper we present Monte Carlo simulations of the 
optically-selected quasar population fine-tuned to study 
the radio loudness distribution of quasars. We investi- 
gate four different models based on single and double 
Gaussian distributions and different radio-to-optical lu- 
minosity relationships. Our aim was to investigate the 
long-standing ambiguity in the existence of the intrin- 
sic difference between radio-quiet and radio-loud quasars 
and in particular, the existence of bimodality in their ra- 
dio loudness distribution. The sample used here, based 
on SDSS DR7 Quasar Catalog matched to the FIRST 
survey, is the largest ever analyzed (8307 radio-detected 
quasars), uniformly selected and reliable (spectroscopi- 
cally confirmed), and the method properly accounts for 
uncertainties in the K-corrections, optical variability and 
survey incompleteness. 

We find that the best-fit model for this SDSS-FIRST 
quasar sample is a two-component model, but the com- 
ponents overlap so that radio loudness bimodality (i.e. a 
minimum between the radio loud and quiet portions of 
the distribution) is not apparent. Statistics of our fits 
indicate that even our best-fit model does not describe 
the data optimally, although it is significantly better than 
any other model we used. Our main result, that the radio 
loudness distribution of quasars likely consists of at least 
two components, agrees with earlier findings indicating 
the existence of two distinct populations of quasars. In 
the framework of the simple two-component models pre- 
sented here, we conclude that bimodality is not likely in a 
sample which, like ours, covers a broad redshift and opti- 
cal luminosity range, even if it is present in narrower red- 
shift and magnitude bins (e.g. at z > 1.5 and i < 18.5). 
We also conclude that the distribution is unlikely to be 
universal, i.e. independent of redshift or optical luminos- 
ity. 

Investigating possible redshift and optical luminosity 
dependence of the radio loudness distribution, we have 
found that a monotonous change of its shape is visible for 
quasars with i < 18.5. It would imply that at high red- 
shift radio-loud quasars were rarer, on average 'louder' 
and less spread out in radio loudness. The same general 
trend is marginally detectable for fainter quasars in the 
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sample, but the smooth change of the radio loudness dis- 
tribution shape is not nearly as apparent as for i < 18.5. 
We further conclude that the radio loudness distribution 
is likely dependent on redshift and/or optical luminosity, 
but we can not disentangle the two dependencies with 
current models and data. We expect more sophisticated 
two-component models to adequately describe the two 
classes of quasars in future work on this topic. 
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